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Abstract 

In this work we propose a Bayesian framework for data fusion of multivariate 
signals which arises in imaging systems. More specifically, we consider the case 
where we have observed two images of the same object through two different 
imaging processes. The objective of this work is then to propose a coherent 
approach to combine these data sets to obtain a segmented image which can 
be considered as the fusion result of these two images. 

The proposed approach is based on a Hidden Markov Modeling (HMM) of 
the images with common segmentation, or equivalently, with common hidden 
classification label variables which is modeled by the Potts Markov Random 
Field. We propose then an appropriate Markov Chain Monte Carlo (MCMC) 
algorithm to implement the method and show some simulation results and ap- 
plications. 
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1 Introduction 

Data fusion and multi-source information has become a very active area of research 
in many domains : industrial nondestructive testing and evaluation ([!]), medical 
imaging [2, 3, 4], industrial inspection ([5]) (quality control and condition monitoring) 
and security systems in general. 

In all these areas, the main problem is how to combine the information contents of 
different sets of multivariate data gi{r). When the data set gi{r) is an image we have 
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r E B?, and the problem becomes how to combine and represent their fusion. Very 
often the data sets g,i do not represent the same quantities. For example in medical 
imaging we have 2D radiographic data gi and echographical data (72 which are related 
to different properties /i and /2 of the body under examination by 

g,{r)^[HMr)+e,{r) (1) 

where Hi are the operator functionnal of the measuring systems. We may note that 
estimating given each set of data gi is an inverse problem by itself which is often 
an ill-posed problem even if we may know perfectly the operator Hi. So very often 
people use the two data sets separately to obtain two images /i and /2 and then they 
try to make a data fusion. We think it is possible to do a better job if we define 
more precisely what we mean by data fusion of two images /i and /2 and if we try to 
use the data gi and ^2 to estimate directly not only /i and /2 but also the common 
feature of them which we present by a third image z. 




Figure 1: Examples of images for data fusion. a,b) Copyright American Science 
and Engineering, Inc., 2003 : two observations from transmission and backscattering 
X rays, c.d) MRI (Magnetic Resonance Imaging) and CT (Computed Tomography) 
images in medical imaging. 

In this paper, to show the same ideas, we consider first the case where the two 
measuring systems can be assumed almost perfect which means that we can write 

9i{r) = fi{r) + ei{r), z = l,2 (2) 

and we focus on defining what we mean by a common feature z of the two images, 
how to model the relation between /j and z and how to estimate /i, /2 and z directly 
from the two data sets gi and g2- 

The applications we have in mind in this work are either medical imaging or security 
systems imaging. As an example of the two data sets in the first application we con- 
sider MRI and CT images and as an example of the second apphcation we consider a 
transmission and a backscattering diffusion images using X rays (see figure 1). 
The rest of the paper is organized as follows : In section 2 we introduce the com- 
mon feature 2;, model the relation between the images fi to it through p{fi\z) and its 
proper characteristics through a prior law p{z), and describe the Bayesian approach 
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to estimate fi, f2 and z through the a posteriori law /2, z\gi,g2)- In section 3 we 
give some details on the selection of a priori probability laws p{0) of the hyperparam- 
eters which define the a posteriori law p(/i, /2, z\gi, g2). In section 4 we give detailed 
expressions of the aforementionned a posteriori law and propose general structure of 
the MCMC algorithm to estimate /i, /2 and z. Finally, in section 5 we present some 
simulation results to show the performances of the proposed method. 

2 Modeling for data fusion 

In this paper we consider the model (2) where after discretization and using the 
notations Qi = [gi{l), gi{S)f, fi =[/,(!),..., fi{S)f and £i = [£^(1), . . . , ei{S)f 
with S the total number of pixels of the images fi, we have : 

gi^fi + £i, i = l,2 (3) 

Within this model and assuming Gaussian independant noises, p{£i) — A/'(0, u^. J), 
we have 

2 2 

p{gi,92\fi, /2) = n^(fl'^i-^^) ^ n^^i^s'^ ~ -^^^ 

1=1 i=l 

As we want to reconstruct an image with statistically homogeneous regions, it is 
natural to introduce a hidden variable z — {z{l) , . . . , z{S)) G {1,..., K}^ which 
represents a common classification of the two images /j. The problem is now to 
estimate the set of variables (/i, f2, z) using the Bayesian approach : 

P{fi,f2,z\gi,g2) = Pifi,f2\z,gi,g2)p{z\gi,g2) 

oc p(</i|/i, z)p{g2\f2, z)p{fi\z)p{f2\z)p{gi\z)p{g2\z)p{z) 

2 

OC piz)Y[p{gi\fi)p{fi\z)p{gi\z) 

i=l 

Thus to be able to give an expression for /2, z\gi,g2) we need to define p{gi\fi), 
P{fi\z), p{gi\z) and p{z). 

Assuming centered, white and Gaussian, we have : 

Pigm = M{f,,alr) 

To assign p{fi\z) we first define the sets of pixels which are in the same class : 

Rk = {r : z{r) = k}, \Rk\ = rik 
fik = {f^ir) : zir) = k} 
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Then we assume that all the pixels of an image fi which are in the same class will be 
characterized by a mean and a variance af^^ : 

PiMr)\z{r) = k) =M{mik,a^k) 
With these notations we have : 



p{fik) = ^{mk'^,<^iki) 

K 



k=l 
K 



The next step is to define p{gi\z). To do this we may use the relation (3) and the 
laws p{fi\z) and p{£i) to obtain 



P{9i{r)\z{r) = k) = Af{mik, a^^^ + aj) 

1 



Finally we have to assign p{z). As we introduced the hidden variable z for finding 
statistically homogeneous regions in images, it is natural to define a spatial depen- 
dance on these labels. The simplest model to account for this desired local spatial 
dependancy is a Potts Markov Random Field model : 



P{z) = exp <^ a ^ 5{z{r) - z{s)) 
^ ' y res seV(r) 

where S is the set of pixels, S{0) = 1, S{t) = si t ^ 0, V(r) denotes the neighbor- 
hood of the pixel r (here we consider a neighborhood of 4 pixels) and a represents 
the degree of the spatial dependance of the variable z. This parameter wxill be stud- 
ied in the next section and fixed for our algorithm. We have now all the necessary 
prior laws p{gi\fi), pifi\z), p{gi\z) and p{z) and then we can give an expression for 
p{fi,f2,z\gi,g2). However these probability laws have in general unknown parame- 
ters such as a^. in p{gi\fi) or rriik and a"^^ in p{fi\z). In a full Bayesian approach, we 
have to assign prior laws to these "hyperparameters". This point is addressed in the 
next section. 



3 Spatial dependance parameter of the Potts model 

In the Potts Markov Random Field (PMRF) model, the parameter a determines the 
spatial dependancy between the pixels. With this model we can expect for controlling 
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the size of the homogeneous regions in images. Indeed, if wc take for example a — 
then we consider the pixels independant, and if we increase the value of a, we increase 
the spatial dependancy. In this section we study with simulations how this spatial 
dependancy involves with the value of a. 



c) a = 0.65 




(d) a = 0.7 




(e) a = 1 



(f) a = 1.3 




Figure 2: Simulations of the PRMF with variation of a 

Figure 2 represents some realisations of the PRMF on (128x128) images of 8 
labels. We can see easily that the size of the homogeneous regions does not increase 
linearly with the value of a. Also the simulations show the existence of a threshold 
between 0.6 and 0.7. Under this threshold the simultations of the PMRF remain 
strongly disturbed. Beyond this threshold some labels become prevalent and the 
homogeneous regions become quickly very large. 

The PMRF model is often used in image processing and gives a possibility to control 
the size of the homogeneous regions. In our algorithm we will fix the value of a to 
1.3 for introducing a strong spatial dependancy between labels. 



4 Prior selection of the hyperparameters 

There is an extensive litterature on the construction of non informative priors. In 
this section we use results from [6] to choose particular priors, taking into account 
the restriction of our particular parmetrical model. In [6] the authors used differential 
geometry tools to construct particular priors. 

Let rrii = {mik)k=i,...,K and cr^ = {o'f^)k=i^,,,^K be the means and the variances of the 
pixels in different regions of the images fi as defined before. We define di as the set 
of all the parameters which must be estimated : 

6>, = (cTj,mi,o-2), i=l,2 
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We choose the prior distribution of 0i of the following form : 

where pe is the likelihood of 0, po is a reference distribution, C is a constant which 
represents the confidence degree we have on pq, Ds is the 6 — divergence and g is the 
Fisher information matrix. 

The authors in ([6]) showed that if we choose this prior distribution for 6i with S = 0, 
we find the conjugate priors. When applied those results for our case, where these 
priors become : 

- Inverse Gamma X^(q;q', /Sq') ^^^1 X^(Q;jo, Ao) respectively for the variances a^, 
and affc, 

- Gaussian J\f{miQ, afg) for the means rriif.. 

The hyper-hyperparameters ai^, /?.jo, rriiQ and a^Q are fixed and the results are not in 
general too sensitive to their exact values. 

5 a posteriori distributions for the Gibbs algo- 
rithm 

The Bayesian approach consists now to estimate the whole set of variables (/i, /2, 2;, ^i, ^2) 
following the joint a posteriori distribution p{fi, f2,z,di,02\gi,g2)- It is difficult 
to simulate a joint sample {fi, f2, ^,61,62) directly from his joint a posteriori dis- 
tribution. However we can note that considering the prior laws defined before, 
we are able to simulate the conditionnal a posteriori laws p{fi, f2, z\gi,g2,di,02) 
and p{di, 92\gi,g2, fi, f2,z). That's why wc propose a Gibbs algorithm to estimate 
{fi: f2, z,6i,62), decomposing this set of variables into two subsets, {fi,f2,z) and 
{Gil 02)- Then the Gibbs algorithm follows : given an initial state (0i,02)^°\ 

Gibbs sampling 

repeat until convergence 

1. simulate (//"'', A^"'', i^"^) ~ /2, 2;|gri, 512, ^^^'^ ^\^2^" 

2. simulate ^Z*"^ ^ p{Oi\gi, ft\ z^'^^) 



We will now define the conditionnal a posteriori distribution we use for the Gibbs 
algorithm. 

sampling 0i\fi,gi,z : 

We have the following relation : 

P{0i\fi,gi,z) (xp{al\fi,gi) p{mi,(Tl\fi,z) 
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and then we can use the Bayes formula : 

Pirui, cr1\fi, z) oc p{fi\z, rrii, a-1)p{mi)p{a-l) 

and 

P{<^l\fi, 9i) oc p{,gi\fi, crl)p{al) 

Those a posteriori distributions are calculated from the prior selection fixed before 
and we have 

- ruiklfi, z, aff., TJiio, ~ J^ii^ik, vf^), with 



'^ik - \ 2 + 2 
^ik ^iO 



'^iklfi^ Ao ~ Wi^ik, Pik), with 
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- (^cilfh 9i ~ Si), with 

Ui — — + aQ, S — total number of pixels 

Sampling /i, /a, z\gi,g2, 0i, 62 : 
Using the Bayes formula we have : 

p{fi, f2, z\gi,g2, Oi, 62) = f2\z,gi,g2, Oi, e2)p{z\gi,g2, Oi, 62) 

Then the sampling of this joint distribution is again obtained through a Gibbs sam- 
phng scheme and then p{fi\gi, z, 6i) by samphng first p{z\gi,g2, di, 62)- For the first 



step we have : 



Piz\gi,g2,Oi,02) oc p{gi,g2\z,0i,O2)piz) 

= p{gi\z,di)p{g2\z,d2)p{z) 



and for the next step we have : 

p{Mrmr),z{r) = k^O,) = N{m,r\^Tr') 
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where 



2apost f Qiir) rriik \ 

As we choosed a Potts Markov Random Field model for the labels, we may note that 
an exact sampling of the a posteriori distribution vip^Q^iQ^^ Gi, O2) is impossible. In 
theory, in each step, we have to implement again a third Gibbs sampling to obtain 
exact samples of z. However this will increase significantly the complexity of the 
algorithm. To obtain a faster algorithm, the solution we propose consists in imple- 
menting only one cycle of the Gibbs sampling for z in each iteration. In fact it comes 
down to decompose the set of variables into three subsets {di, 62), (/i, /2), and z. 
The Gibbs algorithm we propose is then : given an initial state {Oi, 62, 

Gibbs sampling 
repeat until convergence 

1. simulate z^") ~ j9(z|i("-i), £/2, ^V""^\ ^V""^^) 
simulate //""^ ~ p(/i|cjij, z^^^^, 

2. simulate g/""^ ^ //"'', 



As we choosed a first order neighborhood system for the labels, we may also note 
that it is pssible to implement the Gibbs algorithm in parallel. Indeed, we can de- 
compose the whole set of pixels into two subsets forming a chessboard (see figure 2). 
In this case if we fix the black (respectively white) labels, then the white (respectively 
black) labels become independant. 

black labels 



white labels 



Figure 3: Chessboard decomposition of the labels z 

This decomposition reduces the complexity of the Gibbs algorithm because we can 
simulate the whole set of labels in only two steps. The Parallel Gibbs algorithm we 
implemented is then the following : given an initial state {Oi, 62, zY^\ 



2apost 



m. 



apost 
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Parallel Gibbs sampling 



repeat until convergence 

1. simulate is ~ p(2;|zV^''~^\fl'i,fl'2,^V" ^\^V" 
simulate Zw ^"^^ p{z\ Zb , , flfs , ^V" ^\ O2 ) 
simulate ^ p{fi\gi, z^'^\e!''^ 

2. simulate 0/"^ ^ g^Q 

In the following we have implemented this algorithm. 



6 Spatial dependance in the estimated images 

We want now to introduce a dependance between pixels of which are in a same 
homogeneous region. Our aim is to improve the reconstructed images and then (be- 
cause our algorithm is iterative) improve the quality of our classification. We will 
now describe this new modelisation and the modifications it implies. 



6.1 New modelisation on the images fi 

We now consider that pixels /i(r) inside a same region are dependant. However pixels 
being in different regions stay independant. We then introduce a "contour" variable 
q as follows : 

q{r) — if {z{s), s ~ r} are in a same region 
— 1 else 

Then we have the following : 



K 



p{fi\z,q,di) = Ylp{fik\z,q,di) 



k=l 

Let note fif^{r) — {fi{s), s ~ r}. Then we can write : 

p{Mr)\z{r)=k,q{r),fi^{r),di) = Ar(^,, a^) if g(r) = 1 

Note also 

2 

Then we can write the distribution p{gi{r)\z{r) — k, q{r), fiN{f),Oi) '■ 
P{gi{r)\z{r) = k, q{r), fi^ir), 0i) = Af{mf^^r),crl(r) + <^l) 
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6.2 a posteriori distributions 
6.2.1 Sampling f,\gi,z,q,di 

With the same method of section 5, we obtain the a posteriori distribution : 
PiMr)\gi{r),z{r),q{r), fiN{r),6i) = U{maposu (^Ipost), 

with 



2 



9i{r) ruf^^r) 



IT^apost — ^apost I _2 ~'~ ^2 



"Mr) 



2 _ I |_ 



As we choose a spatial dependance between pixels fi{r), we have the same problem 
as for the labels. Indeed an exact sampling of the a posteriori law p{fi\gi,z,q,0i) 
becomes impossible. The solution is, as for the labels, to decompose the set of pixels 
into a chessboard. Let then note /j^y and fi^ rspcctivcly the white and black pixels 
fi{r). Then if we fix /j^, the pixels /^^(r) are independant and we have 

p{fiB\9i,z,q,0i) = Y[ P{fi{r)\9i{r),z{r),q{r),di) 

r black 

and the sy metric relation for fi^^. 

This solution consists then in introducing a Gibbs algortihm for sampling /j and, as 
for the labels, we implement only one cycle of the Gibbs sampling for limiting the 
complexity of our algorithm. 

6.2.2 Sampling z|gi, ^2, /i, /2, g, ^i, ^2 

As for the first case we use the a posteriori distribution : 

Piz\gi,g2, fi, f2, q, 01, O2) oc p{gi\z, fi, q, di)p{g2\z, /s, q, e2)p{z) 

The exact sampling of this distribution is still impossible but we can use the decom- 
position into a chessboard to obtain : 

2 

P{zB\Zw,gNl,gN2: fBl, fB2,QN, 01,02) OC p{zB\zw)Ylp{gNi\z, fiw, Qn, Oi) 

1=1 
2 

= p(2;s|2;iy)J]^ Y[ p{gi{r)\z{r),fNi{r),q{r),Oi) 



i=l r black 



Then we implement for this part one cycle of a Gibbs sampling, as in the first method 
described in section 5. 
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6.2.3 Sampling di\z,gi,fi,q 

We still use the same method to obtain the a posteriori distributions of the parameters 
of 0i. However we have here to decompose the set Rk into to subsets as follows : 

i?fc = -Rfc u rI 

with Rl = {r; z{r) = k,q{r) = i}. Let also note = \Rl\- With this decomposition 
we can calculate the a posteriori distributions of 0i : 



m, 



;k\fi, z, q, ruio, ~ ^fil^^k, ^1k)^ ^ith 



f'ik 




f^j^fe|/i,2;,q,Q;jo,Ao ~2:^(«ife>Afe), with 



, rik 



r^Rl r^Rl 



2 

- (T^. l/j, Qi ~ ig{iyi, Ei), with 

Ui — — + al\ S — total number of pixels 
= l\\9i-fi\\' + Po 
6.3 New Gibbs algorithm 

The only difference between the algorithm of section 5 is the introduction of the new 
variable q and in the sampling of /j. Indeed we have here to decompose fi into two 
subsets and fi^^. The Gibbs algorithm we have implemented is then : 

Parallel Gibbs sampling 

repeat until convergence 

1. simulate i^^") ~ /.[^''^ Q^'^-^), ^2, ^^^'^"'l 4^""'^) 

simulate z^^^^ ~ Q^"-^), ^1, ^2, 0V""'\ ^V""'^) 

2. simulate q*^") ~ 

3. simulate /i"^ ~ ^V""'^) 



simulate fi^ ~ p{fiw\fiB ^^9h 
4. simulate ^/"^ ~ fl^i) 
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7 Simulation and results 



Here we illustrate two applications of the proposed method in cases of medical imaging 
and security systems. The first application is MRI and CT images of a brain which 
are (256 X 256) images. 




Figure 4: Results of data fusion from MRI and CT images. a,b) MRI and CT images 
in medical imaging. c,d) segmentation of the two images taken independantly with 
respectively 6 and 9 labels. e,f) respective reconstruction of the images, g) result of 
data fusion with 9 labels. 

Figure 4 shows the data fusion result of the proposed method comparing independant 
segmentation of the two images and segmentation using data fusion. As it is seen 
on this figure the fusionned segmentation we obtain contains all the regions and 
boundaries of both images. This is particularly visible in the up-center of the image. 
The second application is X ray transmission and backscattering images, which are 
(141 X 192) images. The observed object is a suitcase containing two guns. Figure 5 
shows the result of the proposed method. The independant segmentation of the X ray 
backscttering image show clearly the presence of the right gun, but it is difficult to 
distinguish the left gun whatever the number of labels. In this figure we can see that 
the X ray transmission image brings essential information to precisely distinguish the 
left gun without eliminating the detccton of the one on the right. 
In both applications we have satisfactoring results of image fusion, even when images 
present a great number of homogeneous regions and boundaries. Figures 4 and 5 
show that the proposed method uses both images to increase the performances of the 
segmentation. Note also that the segmentation time of one image independantly or 
as result of the image fusion is practically the same. Indeed the proposed method 
does not really increase the complexity, making fusion and reconstruction in the same 
time. 

However in both cases the reconstructed images are not visibly improved. This is 
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Figure 5: Results of data fusion from X ray images. a,b) two observations from 
transmission and backscattering X rays, c) segmentation of only image a) with 7 
labels, d) segmentation of only image a) with 8 labels. e,f ) respective reconstruction 
of the images, g) result of data fusion with 8 labels. 



mostly due to the assumption that the values of fi{r) at any two different pixels are 

independant. We may expect for better results of reconstruction and segmentation if 
we introduce some local spatial dcpendancy between the neighboring pixels of images 
/i(r). This point is under developpment and we will report soon on the results. 

8 Conclusion 

We proposed a Bayesian method for data fusion of images, with a Potts Markov 
Random Field model on the hidden variable z. We illustrated how the segmentation 
is improved by using data fusion through two applications : MRI and CT images 
in medical imaging and X ray transmission and backscattering images in security 
systems. 

We showed then how reconstruction and fusion can be computed in the same time 
using a MCMC algorithm, which reduce the complexity of the algorithm. 
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